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Abstract. In this paper, we discuss how to efficiently evaluate and assemble general finite 
element variational forms on H(div) and //(curl). The proposed strategy relies on a decomposition 
of the element tensor into a precomputable reference tensor and a mesh-dependent geometry tensor. 
Two key points must then be considered: the appropriate mapping of basis functions from a reference 
element, and the orientation of geometrical entities. To address these issues, we extend here a 
previously presented representation theorem for affincly mapped elements to Piola-mapped elements. 
We also discuss a simple numbering strategy that removes the need to contend with directions of 
facet normals and tangents. The result is an automated, efficient, and easy-to-use implementation 
that allows a user to specify finite element variational forms on H(div) and .ff(curl) in close to 
mathematical notation. 
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1. Introduction. The Sobolev spaces H(div) and H(cm\) play an important 
role in many applications of mixed finite element methods to partial differential equa- 
tions. Examples include second order elliptic partial differential equations, Maxwell's 
equations for electromagnetism, and the linear elasticity equations. Mixed finite ele- 
ment methods may provide advantages over standard H 1 finite element discretizations 
in terms of added robustness, stability, and flexibility. However, implementing iJ(div) 
and -ff (curl) methods requires additional code complexity for constructing basis func- 
tions and evaluating variational forms, which helps to explain their relative scarcity 
in practice. 

The FEniCS project [TSJ [53] comprises a collection of free software components 
for the automated solution of differential equations. One of these components is the 
FEniCS form compiler (FFC) [2Ql ETJ [25]. FFC allows finite element spaces over 
simplicial meshes and multilinear forms to be specified in a form language close to the 
mathematical abstraction and notation. The form compiler generates low-level (CH — h) 
code for efficient form evaluation and assembly based on an efficient tensor contraction. 
Moreover, the FErari project [T9l [22l [23[ [24] has developed specialized techniques for 
further optimizing this code based on underlying discrete structure. FFC relies on 
the Finite element Automatic Tabulator (FIAT) [16j HZl [H] for the tabulation of 
finite clement basis functions. FIAT provides methods for efficient tabulation of finite 
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element basis functions and their derivatives at any particular point. In particular, 
FIAT provides simplicial ff (div) element spaces such as the families of Raviart and 
Thomas [Mj, Brezzi, Douglas, and Marini [TU], and Brezzi et al. [5], as well as if (curl) 
elements of the Nedelec types [29, 30 . 

Previous iterations of FFC have enabled easy use of ff 1 and L 2 conforming fi- 
nite element spaces, including discontinuous Galerkin formulations, but support for 
_ff(div) and ff(curl) spaces has been absent. In this paper, we extend the previous 
work [201 [HI [31] to allow simple and efficient compilation of variational forms on 
ii(div) and ff(curl), including mixed formulations on combinations of ff , ff(div), 
ff(curl), and L 2 . The efficiency of the proposed approach relies, in part, on the tensor 
representation framework established in [2T]. In this framework, the element tensor 
is represented as the contraction of a reference tensor and a geometry tensor. The 
former can be efficiently precomputed given automated tabulation of finite element 
basis functions on a reference element, while the latter depends on the geometry of 
each physical element. For this strategy, a key aspect of the assembly of ff(div) 
and ff (curl) conforming element spaces becomes the Piola transformations, isomor- 
phically mapping basis functions from a reference element to each physical element. 
Also, the orientation of geometrical entities such as facet tangents and normals must 
be carefully considered. 

Implementations of ff(div) and ff(curl) finite element spaces, in particular of 
arbitrary degree, are not prevalent. There are, to our knowledge, no implementations 
that utilize the compiled approach to combine the efficiency of low-level optimized 
code with a fully automated high-level interface. Some finite element packages, such 
as FEAP [T], do not provide ff(div) or ff(curl) type elements at all. Others, such 
as FreeFEM [35], typically provide only low-order elements such as the lowest-order 
Raviart-Thomas elements. Some libraries such as deal. II [8] or FEMSTER [12] do 
provide arbitrary degree elements of Raviart-Thomas and Nedelec type, but do not 
automate the evaluation of variational forms. NGSolve |36] provides arbitrary order 
ff (div) and -ff (curl) elements along with automated assembly, but only for a prede- 
fined set of bilinear forms. 

This exposition and the FFC implementation consider the assembly of ff (div) and 
ff(curl) finite element spaces on simplicial meshes. However, the underlying strategy 
is extendible to nonsimplicial meshes and tensor-product finite element spaces defined 
on such meshes. A starting point for an extension to isoparametric ff 1 conforming 
finite elements was discussed in [30]. The further extensions to ff(div) and ff(curl) 
follow the same lines as for the simplicial case discussed in this note. 

The outline of this paper is as follows. We begin by reviewing basic aspects of the 
function spaces ff(div) and ff(curl) in section [21 and we provide examples of varia- 
tional forms defined on these spaces. We continue, in section [3j by summarizing the 
ff(div) and ff(curl) conforming finite elements implemented by FIAT. In section 2] 
we recap the multilinear form framework of FFC and we present an extension of the 
representation theorem from [21] . Subsequently, in section [5j we provide some notes 
on the assembly of ff(div) and ff(curl) elements. Particular emphasis is placed on 
aspects not easily found in the standard literature, such as the choice of orientation 
of geometric entities. In section [6] we return to the examples introduced in section [2] 
and illustrate the ease and terseness with which even complicated mixed finite element 
formulations may be expressed in the FFC form language. Convergence rates in agree- 
ment with theoretically predicted results are presented to substantiate the veracity of 
the implementation. Finally, we make some concluding remarks in section [7] 
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2. i?(div) and iJ(curl). In this section, we summarize some basic facts about 
the Sobolev spaces if(div) and -ff(curl) and we discuss conforming finite element 
spaces associated with them. Our primary focus is on properties relating to interele- 
ment continuity and change of variables. The reader is referred to the monographs [11) 
and [28] for a more thorough analysis of H(div) and If (curl), respectively. 

2.1. Definitions. For an open domain SI C R™, we let L 2 (S1,R") denote the 
space of square- integrable vector fields on SI with the associated norm 1 1 • | |o and inner- 
product (•,•), and we abbreviate L 2 (S1) = L 2 (Q, R 1 ). We define the following standard 
differential operators on smooth fields v: D a v = S? 1 • • • d" m v for a multi-index a of 
length m, divv = Yn=i dxi v i, m ^ v = {d X2 v z - d X3 v 2 ,d X3 vi - d Xl v 3l d Xl v 2 - d X2 v{), 
and rotu = d Xl v 2 — d X2 v\. We may then define the spaces H" m (Sl), iJ(div;Sl), and 
ff(curl; f2) by 

{v G L 2 (S1) : D a v G L 2 (fl), \a\ < m}, m = 1, 2, . . . , 

{v G L 2 (n,« n ) : divv G L 2 (n)}, 

f {v G L 2 (n,M. 2 ) : rotv G L 2 (fl)}, !!cM 2 , 
\ {v G L 2 (17,R 3 ) : curlw G L 2 (f2,R 3 )}, !!cM 3 , 

with derivatives taken in the distributional sense. The reference to the domain Q 
will be omitted when appropriate, and the associated norms will be denoted || • || m , 
| • ||div, and || • 1 1 C uri- Furthermore, we let M denote the space of matrices and 
we let if(div; f2,M) denote the space of square-integrable matrix fields with square- 
integrable row- wise divergence. 

For the sake of compact notation, we shall also adopt the exterior calculus nota- 
tion of [5] and let A fc (f2) denote the space of smooth differential /c-forms on SI, and let 
L 2 A k (tt) denote the space of square-integrable differential k- forms on fl. We further 
let d denote the exterior derivative with adjoint 6, and we define HA k (il) = {v G 
L 2 A k (il), dti 6 L 2 A k (il)}. Further, V r A k is the space of polynomial fc-forms of up to 
and including degree r, and V~A k denotes the reduced space as defined in section 3.3] . 

2.2. Examples. The function spaces iJ(div) and -ff(curl) are the natural func- 
tion spaces for an extensive range of partial differential equations, in particular in 
mixed formulations. We sketch some examples in the following, both for motivational 
purposes and for later reference. The examples considered here are mixed formulations 
of the Hodge Laplace equations, the standard eigenvalue problem for Maxwell's equa- 
tions, and a mixed formulation for linear elasticity with weakly imposed symmetry. 
We return to these examples in section [6] 

Example 2.1 (mixed formulation of Poisson's equation). The most immediate 
example involving the space if (div) is a mixed formulation of Poisson's equation: 
—Am = / in fl C W 1 . By introducing the flux a = — gradw, and assuming Dirichlet 
boundary conditions for u, we obtain the following mixed variational problem: Find 
a G H(div; f2) and u G L 2 (J1) satisfying 

(2.1) (r, a) - (div r, u) + (v, div a) = (v, f) 

for all r G #(div; fi) and v G L 2 (n). 

Example 2.2 (the Hodge Laplacian) . With more generality, we may consider weak 
formulations of the Hodge Laplacian equation (d(5 + <5d)w = / on a domain fl C R n ; 
see [51 section 7] . For simplicity of presentation, we assume that SI is contractible such 
that the space of harmonic forms on SI vanishes. The formulation in Example 12.11 is 



H(div; f2) = 
F(curl; Q) = 



ASSEMBLY OF H(div) AND tf(curl) FINITE ELEMENTS 



4133 



the equivalent ol seeking u G HA n and a = Su G i/A™ -1 for n = 2,3 with natural 
boundary conditions (the appropriate trace being zero). To see this, we test a — Su 
against r G HA n ^ 1 and we test (d 8 + S d)u = / against v 6 HA n to obtain 

(r,er) - (t,5u) + (v,da) = (v,f), 

noting that du — for u G HA n . Integrating by parts, we obtain 

(2.2) (T,a)-(dT,u) + (v,da) = («,/>. 

We may restate (|2.2I) in the form (|2.ip by making the identifications <5u = — grad it, 
dr = divr, and der = diver. If Q C R 3 , we may also consider the following mixed 
formulations of the Hodge Laplace equation. 

(i) Find a G HA 1 = if (curl) and u G HA 2 = iJ(div) such that 

(2.3) (r, a) — (curlr, u) + (v, curler) + (divw, divu) = (v,f) 

for all r G HA 1 , v G HA 2 . 

(ii) Find a G H A = F 1 and u G i/A 1 = iJ(curl) such that 

(2.4) (r, a) — (grad r, u) + (v, grad tr) + (curl v, curl u) = («, /) 

for all r G i/A°,w G -ffA 1 . 
Example 2.3 (cavity resonator). The time-harmonic Maxwell equations in a 
cavity with perfectly conducting boundary induces the following eigenvalue problem: 
Find resonances uj G K and eigenfunctions E G 7?o(curl; ft), satisfying 

(2.5) (curl F, curl E) = u 2 (F,E) VF G H (cml;Cl), 

where i?o(curl; ft) = {v G i/(curl; 0)| u X = 0}. Note that the formulation (|2.5I) 
disregards the original divergence-free constraint for the electric field E and thus 
includes the entire kernel of the curl operator, corresponding to lj = and electric 
fields of the form E — grad ip. 

Example 2.4 (elasticity with weakly imposed symmetry). Navier's equations for 
linear elasticity can be reformulated using the stress tensor <r, the displacement u, 
and an additional Lagrange multiplier 7 corresponding to the symmetry of the stress 
constraint. The weak equations for il C R 2 , with the naturajj boundary condition 
u\ m = 0, take the following form: Given / G L 2 (f7,R"), find a G iJ(div;0,M), 
u G L 2 {fl,W l ), and 7 G L 2 (tt) such that 

(2.6) (r, Act) + (div r, u) + (v, div a) + (skw r, 7) + (77, skw a) = (v, f) 

for all t G ff(div;n,M), u G L 2 (r2,R"), and 7? G i 2 (0). Here, A is the compliance 
tensor, and skwr is the scalar representation of the skew-symmetric component of r; 
more precisely, 2skwr = ti\— t\%. This formulation has the advantage of being robust 
with regard to nearly incompressible materials and it provides an alternative founda- 
tion for complex materials with nonlocal stress-strain relations. For more details, we 
refer the reader to [B]. 



1 Note that the natural boundary condition in this mixed formulation is a Dirichlet condi- 
tion, whereas for standard H 1 formulations the natural boundary condition would be a Neumann 
condition. 
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2.3. Continuity-preserving mappings for iJ(div) and .ff(curl). At this 
point, we turn our attention to a few results on continuity-preserving mappings for 
H(div) and -ff(curl). The results are classical and we refer the reader to [TTJ [55] for 
a more thorough treatment. 

First, it follows from Stokes' theorem that in order for piecewise iJ(div) vector 
fields to be in H(drv) globally, the traces of the normal components over patch in- 
terfaces must be continuous, and analogously tangential continuity is required for 
piecewise -ff(curl) fields. More precisely, we have the following: Let Th — {K} be a 
partition of f2 into subdomains. Define the space E/j of piecewise H (div) functions 
relative to this partition Th- 

(2.7) E,, L 2 (Q,R n ) : <j>\ K G H (div; K) G %}■ 

Then <f> G E^ is in H (div; fi) if and only if the normal traces of cj> are continuous 
across all element interfaces. Analogously, if 4>\k G iJ(curl; K) for all K £ Th, then 
4> G iJ(curl; ft) if and only if the tangential traces are continuous across all element 
interfaces. 

Second, we turn to consider a nondegenerate mapping F : Qq — > F(Qq) — Q with 
Jacobian DF(X), X G tt G R". For $ G H m (n Q ), the mapping T dehned by 

(2.8) 7($) = $oF" 1 

is an isomorphism from H m (ilo) to H m (fl). This, however, is not the case for £T(div) 
or _ff(curl), since T does not in general preserve continuity of normal or tangential 
traces. Instead, one must consider the contravariant and covariant Piola mappings 
which preserve normal and tangential continuity, respectively. 

Definition 2.5 (the contravariant and covariant Piola mappings). Let f2 C R n , 
let F be a nondegenerate mapping from Q onto F(Qq) — with J — DF(X), and 
let $ G L 2 (ft ,R"). 

The contravariant Piola mapping J-" dlv is defined by 

(2.9) J- div ($) = — J^oT 1 . 

det J 

The covariant Piola mapping J rcurl i$ defined by 

(2.10) j tmI ($) = r T $or 1 . 

Remark 2.6. We remark that the contravariant Piola mapping is usually dehned 
with an absolute value, J rdlv (<i>) = . dc 1 t j, Jj> o F^ 1 . However, omitting the absolute 

value, as in (|2.9[) . can simplify the assembly of H(div) elements, as will be expounded 
in section [5j 

The contravariant Piola mapping is an isomorphism of H (div; fio) onto iJ(div; Q), 
and the covariant Piola mapping is an isomorphism of -ff (curl; 51o) onto iJ(curl; tt). In 
particular, the contravariant Piola mapping preserves normal traces, and the covari- 
ant Piola mapping preserves tangential traces. We illustrate this below in the case 
of simplicial meshes in two and three space dimensions (triangles and tetrahedra). 
The same results hold for nonsimplicial meshes with cell-varying Jacobians, such as 
quadrilateral meshes [H [TT] . 

Example 2.7 (Piola mapping on triangles in R 2 ). Let Ko be a triangle with 
vertices X 1 and edges E l for i = 1, 2, 3. We define the unit tangents by T l = E t /\\E l \\. 
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We further define the unit normals by N % = RT l , where 

(2.H) H=(_° 1 J 

is the clockwise rotation matrix. 

Now, assume that K Q is affinely mapped to a (nondegenerate) simplex K with 
vertices x 1 . The afhne mapping Fk ■ K — > K takes the form x = Fk(X) = JX + b 
and satisfies x l — F K (X l ) for i = 1, 2, 3. It follows that edges are mapped by 

e = x l - = JpT - X J ) = XE. 

Similarly, normals are mapped by 

||e||n = Re = RJE = (det J)J~ T RE = (det J) J- T \\E\\N, 

dcE7 J 



where we have used that -^t-jRJR t — J T and thus RJ = (det J) J T R for J E 



p2x2 



The relation between the mappings of tangents and normals (or edges and rotated 
edges) may be summarized in the following commuting diagrams: 



(2.12) 



T JWEW/We^ t E _/_^ e 

r\ r R R 



J-T 



N (dct.7).7- r ||^||/||e|| ) n mN (dot .7)7- > 



With this in mind, we may study the effect of the Piola transforms on normal 
and tangential traces. Let $ G C°°(K a , W l ) and let = 7 rdiv ($). Then 

||e|| <j>(x) - = 1 1 e 1 1 ((dot J)~ 1 J<i>(X)) T ((det J) ,r T \ \E\ \/\ \e\ \N) = \ \E\\ • N. 

Thus, the contravariant Piola mapping preserves normal traces for vector fields under 
affine mappings, up to edge lengths. In general, the same result holds for smooth, 
nondegenerate mappings Fk if the Jacobian DFk(X) is invertible for all X G Kq. 
Similarly, let 4> = F CU1 \<S>). Then 

(2.13) \\eU(x) -t=\\e\\ (J- T <f(X)f (J||£||/|M|) = \\E\\ $(X) • T. 

Thus, the covariant Piola preserves tangential traces for vector fields, again up to edge 
lengths. Observe that the same result holds for tetrahedra without any modifications. 
The effect of the contravariant and covariant Piola mappings on normal and tangential 
traces is illustrated in Figure I2TT1 where \ \E\ \ = ||e|| for simplicity. 

Example 2.8 (contravariant Piola mapping on tetrahedra in R 3 ). Now, let Kq be 
a tetrahedron. As explained above, the covariant Piola mapping preserves tangential 
traces. To study the effect of the contravariant Piola mapping on normal traces, we 
define the face normals of K by N — ||^* . Then 

He' x e j \\n = JE l x JE j = det JJ~ T {E l x E 3 ) = \\E l x E-' \ \ det JJ~ T N, 

since (Ju) x (Jv) = det JJ~ T {u x v). Let $ G C°°(K Q ,W l ) and let = J" div ($). 
Then, it follows that 

\\e l x e j \\(t>{x) - n= \\E l x E 3 \\$(X) ■ N. 
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Fig. 2.1. Mapping two vector fields <3> n and <3?t between two triangles using the contravariant 
and covariant Piola mappings. The contravariant Piola mapping (above) preserves normal traces 
of vector fields, and the covariant Piola mapping (below) preserves tangential traces of vector fields. 
This means in particular that the contravariant Piola mapping maps tangents to tangents (which 
have a zero normal component), and that the covariant Piola mapping maps normals to normals 
(which have a zero tangential component). Note that this is somewhat counterintuitive; the con- 
travariant H(div) Piola mapping always maps tangential fields to tangential fields but does not in 
general map normal fields to normal fields. However, in both cases the normal component (being 
zero and one, respectively) is preserved. 

Thus, the contravariant Piola mapping preserves normal traces, up to the area of 
faces. 

We finally remark that if J G R 2x2 defines a conformal, orientation-preserving 
map, the contravariant and covariant Piola mappings coincide. In R 3 , J must also be 
orthogonal for this to occur. 

3. H(div) and i?(curl) conforming finite elements. To construct -ff(div) 
and iJ(curl) conforming finite element spaces, that is, discrete spaces Vh satisfying 
Vh C iJ(div) or Vh C ff(curl), one may patch together local function spaces (finite 
elements) and make an appropriate matching of degrees of freedom over shared ele- 
ment facets. Here, a facet denominates any geometric entity of positive codimension 
in the mesh (such as an edge of a triangle or an edge or face of a tetrahedron). In 
particular, one requires that degrees of freedom corresponding to normal traces match 
for ff(div) conforming discretizations and that tangential traces match for _ff(curl) 
conforming discretizations. 

Several families of finite element spaces with degrees of freedom chosen to facilitate 



ASSEMBLY OF H(div) AND tf(curl) FINITE ELEMENTS 



4137 



this exist. For H(div) on simplicial tessellations in two dimensions, the classical 
conforming families are those of Raviart and Thomas (RT r , r = 0,1,2,...) [3~i] ; 
Brezzi, Douglas, and Marini (BDM r , r = 1,2, ...) [TDJ; and Brezzi et al. (BDFM,., 
?' = 1,2,...) [9]. The former two families were extended to three dimensions by 
Nedelec [29] [30]. However, the same notation will be used for the two- and three- 
dimensional H (div) element spaces here. For iJ(curl), there are the families of Nedelec 
of the first kind (NED r , r = 0, 1, 2, . . . ) [29] and of the second kind (NED;;, r = 
1,2,...) [3D]. We summarize in Table I3TT1 those H (div) and i/(curl) conforming finite 
elements that are supported by FIAT and hence by FFC. In general, FFC can wield 
any finite element space that may be generated from a local basis through either of the 
aforedescribed mappings. In Table I3T21 we also summarize some basic approximation 
properties of these elements for later comparison with numerical results in section [BJ 

Table 3.1 

H(div) and H(cuvl) conforming finite elements on triangles and tetrahedra supported by FIAT 
and FFC for r > 1. When applicable, the elements are listed with their exterior calculus no- 
tation, along with their original references. Note that for K C K 3 , the Raviart-Thomas and 
Brezzi-Douglas-Marini elements are also known as the first and second kind H(div) Nedelec (face) 
elements, respectively. 



Simplex 


iT(div) 


_ff(curl) 


K CM 2 


BDM r VrA 1 ^) [10] 

RT r _i Tr^{K) [34] 
BDFM r [9] 


NED r _i 


ffCK 3 


BDM r V r A- 2 (K) [30] 
RT r _i Vr^ 2 (K) [29] 
BDFM r 


NED r _ 1 VrA^K) [29] 



Table 3.2 

Approximation properties of the spaces from Table [3~TT C > 0, r > X, 1 < m < r. Tl^ denotes 
the canonical interpolation operator, defined by the degrees of freedom, onto the element space in 
question. For simplicity, it is assumed that v is sufficiently smooth for the interpolation operators 
to be well-defined, and for the given norms to be bounded. For more details and sharper estimates; 

cf. unm]. 



Finite element 


Interpolation estimates 


VrA k (Q.) 


\\V-U h v\\ < Ch^Mm + u 


\\v - n h v\ 


Idiv.curl < Cfc m ||j;|| m +1 




\\v-U h v\\ < Ch m \\v\\ m , 


\\v - U h v\ 


div, curl < Ch m l)|| m +l 


BDFMr 


\\v-n h v\\ < ch m \\v\\ m , 


\\v - n h v\ 


Idiv < Ch™\\v\\ m + 1 



For the reasons above, it is common to define the degrees of freedom for each 
of the elements in Table 13.11 as moments of either normal or tangential traces over 
element facets. However, one may alternatively consider point values of traces at 
suitable points on element facets (in addition to any internal degrees of freedom). 
Thus, the degrees of freedom for the lowest order Raviart-Thomas space on a triangle 
may be chosen as the normal components at the edge midpoints, and for the lowest 
order Brezzi-Douglas-Marini space, we may consider the normal components at two 
points on each edge (positioned symmetrically on each edge and not touching the 
vertices). This, along with the appropriate scaling by edge length, is how the degrees 
of freedom are implemented in FIAT. 
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4. Representation of if (div) and Ji(curl) variational forms. In this sec- 
tion, we discuss how multilinear forms on H (div) or H (curl) may be represented as 
a particular tensor contraction, allowing for precomputation of integrals on a refer- 
ence element and thus efficient assembly of linear systems. We follow the notation 
from [20] [21] and extend the representation theorem from [21] for multilinear forms 
on H 1 and L 2 to H(div) and -ff(curl). The main new component is that we must use 
the appropriate Piola mapping to map basis functions from the reference element. 

4.1. Multilinear forms and their representation. Let f2 C M™ and let 

{V^}j =1 be a set of finite dimensional spaces associated with a tessellation T = {K} 
of fl. We consider the following canonical linear variational problem: Find uu € V 2 
such that 

(4.1) a(v,u h )=L(v) VveV,l 

where a and L are bilinear and linear forms on x Vfi and V£, respectively. Dis- 
cretizing (|4.ip . one obtains a linear system AU = b for the degrees of freedom U of 
the discrete solution Uh- 

In general, we shall be concerned with the discretization of a general multilinear 
form of arity p, 

(4.2) a : V£ x V 2 x • • • x V£ -> R. 

Typically, the arity is p = 1 (linear forms) or p = 2 (bilinear forms), but forms of higher 
arity also appear (see 20 ). For illustration purposes, we consider the discretization 
of the mixed Poisson problem (|2.1[) in the following example. 

Example 4.1 (discrete mixed Poisson). Let and Wh be discrete spaces approx- 
imating H(div; f2) and L 2 (Q), respectively. We may then write (|2.1[) in the canonical 
form (|4.ip by defining 

(4.3a) a((r h ,v h ), (a h ,u h )) = (r h ,(7 h ) - (divT h ,u h ) + (v h ,diva h ), 

(4.3b) L((r h ,v h )) = (v h J) 

for (r h ,v h ) e V,l = S h x W h and (a h ,u h ) € V£ = V%. 

To discretize the multilinear form (|4.2[) . we let {(f> 3 k }^l 1 denote a basis for Vl for 
j = 1, 2, . . . , p and we define the global tensor 

(4-4) A- $,,...,#,), 

where i = . . . ,ip) is a multi-index. Throughout, j, k denote simple indices. If 

the multilinear form is defined as an integral over D, — UxeTh^t the tensor A may be 
computed by assembling the contributions from all elements, 

(4-5) At =a((pj 1 ,4> 2 2 ,...,(j)^) = J2 a K (^,^ 2 ,...,^), 

Ken 

where a K denotes the contribution from element K. We further let {(/) k ' J }2Li denote 
the local finite element basis for V£ on K and define the element tensor A K by 

(4-6) Af = a K {^\^\...,<t>?;»). 
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The assembly of the global tensor A thus reduces to the computation of the element 
tensor A K on each element K and the insertion of the entries of A K into the global 
tensor A. 

In [21] , it was shown that if the local basis on each element K may be obtained as 
the image of a basis on a reference element Kq by the standard (affine) isomorphism 
Tk '■ FJ y (Kq) — >• i? 1 (i4T), then the element tensor A K may be represented as a tensor 
contraction of a reference tensor A , depending only on the form a and the refer- 
ence basis, and a geometry tensor Gk, depending on the geometry of the particular 
element K, 

(4-7) Af — A® a G K , 

with summation over the multi-index a. It was further demonstrated in |21j that 
this representation may significantly reduce the operation count for computing the 
element tensor compared to standard evaluation schemes based on quadrature. 

Below, we extend the representation (I4.7[) to hold not only for bases that may be 
affincly mapped from a reference element, but also for finite element spaces that must 
be transformed by a Piola mapping. 

4.2. A representation theorem. We now state the general representation the- 
orem for multilinear forms on H 1 , H(cml), H(div) (and L 2 ). Instead of working out 
the details of the proof here, we refer the reader to the proof presented in [2T] for H 1 , 
and we illustrate the main points for H(div) and -ff (curl) by a series of examples. 

Theorem 4.2. Let Kq C R" be a reference element and let F K : Kq — > K = 
Fk(Ko) be a nondegenerate, affine mapping with Jacobian Jk- For j = 1,2,...,/}, let 
{<f> k denote a basis on K generated from a reference basis {$t.}fc on Kq, that is, 
<jy£'3 = ^Ffci^k)' where F J K is either of the mappings defined by (|2.8[) . (|2.9[) . or (|2.10[) . 

Then there exists a reference tensor A® , independent of K , and a geometry tensor 
Gk such that A K = A : Gk, that is, 

(4.8) Af = Y J A° ta G K Wei, 

for a set of primary indices T and secondary indices A. In fact, the reference tensor 
A takes the following canonical form: 

(4-9) 4 a = Y,i n^^IOldX; 

that is, it is the sum of integrals of products of basis function components and their 
derivatives on the reference element Kq, and the geometry tensor Gk is the outer 
product of the coefficients cr.) of any weight functions with a tensor that depends only 
on the Jacobian Jk, 

(aim r a -TJr |det Jk| STT[ dX U TT 9x U 

(4.10) Gk _ 11c( . ) __^ 11 __ 11 __, 

for some integer 7. 

4.3. Examples. To this end, we start by considering the vector- valued L 2 (il) 
inner product, defining a bilinear form: 

(4.11) a(v, u) = / v-udx. 

Jn 



4140 



MARIE E. ROGNES, ROBERT C. KIRBY, AND ANDERS LOGG 



In the following, we let x denote coordinates on K and we let X denote coordi- 
nates on the reference element K . Fx is an affine mapping from K a to K, that is, 
x = Fk{X) = JkX + xk- We further let <j> K denote a held on K obtained as the 
image of a field $ on the reference element Kq, (f> K = jty (<$>). We aim to illustrate 
the differences and similarities of the representations of the mass matrix for different 
choices of mappings J-k, in particular, affine, contravariant Piola, and covariant Piola. 

Example 4.3 (the mass matrix with afhnely mapped basis). Let Tk be the affine 
mapping, Fk{Q) = $ o F% . Then, the element matrix A K for (I4.11[) is given by 

(4.12) Af=f ^ x {x)-<j>^{x)dx=\detJ K \ [ ^ \ff\(X) ^M(X) dX, 

JK Jk 

where we use $[/3] to denote component /3 of the vector- valued function $ and implicit 
summation over the index (3. We may thus represent the element matrix as the tensor 
contraction (|4.7[) with reference and geometry tensors given by 

A° i= [ <S>l[p](X)S>l[0](X)dX, 
Jk 

G K = |det J K \. 

We proceed to examine the representation of the mass matrix when the basis 
functions are transformed with the contravariant and the covariant Piola transforms. 

Example 4.4 (the mass matrix with contravariantly mapped basis). Let J~x v be 
the contravariant Piola mapping, 

^ V($) = de^ J ^ OF - 
Then, the element matrix A K for (|4.1ip is given by 

Af= [ ^ 1 {x).^\x)dx 



\det Jk\ dxp dxp 



(det J K ) 2 dX ai dX, 



a 2 J K 



$l[ ai }(X)<S>l[a 2 }(X)dX. 



We may thus represent the element matrix as the tensor contraction (|4.7I) with refer- 
ence and geometry tensors given by 



Al = 



[ ^ i [a 1 ](X)$2 [a2 ]( X )dX, 
Jk 



K _ |detJjf| dxp dxfj 
a ~ (det J K f 3X ai dX a2 ' 

Example 4.5 (the mass matrix with covariantly mapped basis). Let T^ 11 be the 
covariant Piola mapping, 

JT 1 ($) = J/$oF k 1 . 
Then, the element tensor (matrix) A K for (|4. 1 1 [) is given by 

Af= j ^\x)-^ 2 (x)dx 

JK 

-l dEt *ll7^?l*- |£,ll(x) ' t » H(A ' )dA '- 
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We may thus represent the element matrix as the tensor contraction (|4.7I) with refer- 
ence and geometry tensors given by 



A° ia = I ^UaiKX) 3> 2 i2 [a 2 }(X)dX, 

JKn 



G* = I det J K 



K„ 

dX a , dX a 



We observe that the representation of the mass matrix differs for affine, con- 
travariant Piola, and covariant Piola. In particular, the geometry tensor is different 
for each mapping, and the reference tensor has rank two for the affine mapping, 
but rank four for the Piola mappings. We also note that the reference tensor for 
the mass matrix in the case of the covariant Piola mapping transforms in the same 
way as the reference tensor for the stiffness matrix in the case of an affine mapping 
(see |H]). 

It is important to consider the storage requirements for this tensor contraction 
approach and when other approaches might be appropriate. For either the -ff(div) 
or iJ(curl) mass matrix, for example, the reference tensor A has rank four (two 
indices for vector components and two for basis functions). As such, the storage 
requirements for A are cPn 2 , where d = 2, 3 is the spatial dimension and n is the 
number of reference element basis functions. We also note that n = 0(r d ), where r 
is the polynomial degree. Storing A is thus comparable to storing d 2 element mass 
matrices. This is a modest, fixed amount of storage, independent of the mesh. The 
tensor contraction may be computed in several different ways. The default option 
used by FFC is to generate straightline code for performing the contraction of A 
and G . Alternatively, one may also consider A being stored in memory as an array 
and applied via BLAS. In the first case, the size of generated code can become a 
problem for complex forms or high-order methods, although this is not as large of a 
problem in the second case. The geometry tensor, G , must be computed for each 
element of the mesh. For either the contravariant or covariant case, G K is a d x d 
array and so is comparable to storing the cell Jacobian for each cell of the mesh. 
For more complicated forms, storing G K for each cell can become more expensive. 
However, FFC currently stores only one such G K at a time, interleaving construction 
of G K and its multiplication by A . For more complicated bilinear forms (such as 
ones involving multiple material coefficients), the memory requirements of A and 
G K both grow with the polynomial degree, which can lead to inefficiency relative to a 
more traditional, quadrature-based approach. For a thorough study addressing some 
of these issues, we refer the reader to [32] . 

FFC is typically used to form a global sparse matrix, but for high-degree ele- 
ments, static condensation or matrix-free approaches will be more appropriate. This 
is a result of the large number of internal degrees of freedom being stored in the 
sparse matrix and is an artifact of assembling a global matrix rather than our tensor 
contraction formulation as such. 

We conclude by demonstrating how the divergence term from (|4.3[) is transformed 
with the contravariant Piola (being the relevant mapping for if(div)). 

Example 4.6 (divergence term). Let Tk be the affine mapping, let be the 
contravariant Piola mapping, and consider the bilinear form 

(4.13) a(v 7 a) = / vdivcrdx 

Jk 
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for (v,a) e V 1 x V 2 . Then, if <j) K ' 1 = J r K ($ 1 ) and (f> K - 2 = Jf v (<I> 2 ), the element 
matrix A K for (|4.13[) is given by 

7x 1 2 detJ K <9A Ql ftr^ J Ko 11 dX a2 

Noting that q^ 13 d gJ" 2 — S aia21 we may simplify to obtain 

We may thus represent the element matrix as the tensor contraction (|4.7I) with refer- 
ence and geometry tensors given by 

A° = / div$ 2 2 dA, 
G K a = ±1. 

The simplification in the final example is a result of the isomorphism, induced 
by the contravariant Piola transform, between i?(div, Kq) and H (div, K). FFC takes 
special care of such and similar simplifications. 

5. Assembling iJ(div) and i?(curl) elements. To guarantee global conti- 
nuity with Piola-mapped elements, special care has to be taken with regard to the 
numbering and orientation of geometric entities, in particular the interplay between 
local and global orientation. This is well known, but is rarely discussed in the stan- 
dard references, though some details may be found in [23|35]. We discuss here some 
of these issues and give a strategy for dealing with directions of normals and tan- 
gents that simplifies assembly over ff(div) and -ff(curl). In fact, we demonstrate 
that one may completely remove the need for contending with directions by using an 
appropriate numbering scheme for the simplicial mesh. 

5.1. Numbering scheme. The numbering and orientation of geometric entities 
in FFC follows the UFC specification [3J. In short, the numbering scheme works as 
follows. A global index is assigned to each vertex of the tessellation Th (consisting 
of triangles or tetrahedra). If an edge adjoins two vertices Vi and Vj, we define the 
direction of the edge as going from vertex Vi to vertex Vj if i < j. This gives a unique 
orientation of each edge. The same convention is used locally to define the directions 
of the local edges on each element. Thus, if an edge adjoins the first and second 
vertices of a tetrahedron, then the direction is from the first to the second vertex. A 
similar numbering strategy is employed for faces. The key now is to require that the 
vertices of each element are always ordered based on the their global indices. 

For illustration, consider first the two-dimensional case. Let Kq be the UFC 
reference triangle, that is, the triangle defined by the vertices {(0, 0), (1, 0), (0, 1)}. 
Assume that K = Fk{Kq) and K' — Fk>(Kq) are two physical triangles sharing an 
edge e with normal n. If e adjoins vertices Vi and Vj and is directed from Vi to Vj, it 
follows from the numbering scheme that i < j. Since the vertices of both K and K' 
are ordered based on their global indices, and the local direction (as seen from K or 
K') of an edge is based on the local indices of the vertices adjoining that edge, this 
means that the local direction of the edge e will agree with the global direction, both 
for K and K' . Furthermore, if we define edge normals as clockwise rotated tangents, 
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Fig. 5.1. Two adjacent triangles will always agree on the direction of a common edge tangent 
or normal. The two triangles in the figure share a common edge between the global vertices V20 
and v$o- These two vertices have different local indices for K (1,2) and K 1 (2,3), but the ordering 
convention, local numbering according to ascending global indices, ensures that both triangles agree 
on the direction of the common edge e. 



K and K' will agree on the direction of the normal of the common edge. The reader 
is encouraged to consult Figure ISTTI for an illustration. 

The same argument holds for the direction of edges and face normals in three 
dimensions. In particular, if face normals are consistently defined in terms of edges, 
it is straightforward to ensure a common direction. Consider two tetrahedra K and 
K' sharing a face /, defined by three vertices , u,- 2 , i>j 3 such that i\ < 12 < H- 
Clearly, Vi 1 will be the vertex with the lowest index of the face / for both K and K' . 
Furthermore, each of the two edges that adjoin , that is, the edge from Vi t to vi 2 
and the edge from to Vi 3 , has a unique direction by the previous arguments. These 
two edges can therefore define consistent tangential directions of the face. Taking the 
normalized cross-product of these edges gives a consistent face normal. This is the 
approach used by FIAT/FFC. As a consequence, two adjacent tetrahedra sharing a 
common face will always agree on the direction of the tangential and normal directions 
of that face. This is illustrated in Figure [5721 

We emphasize that the numbering scheme above does not result in a consistent 
orientation of the boundary of each element. It does, however, ensure that two adja- 
cent elements sharing a common edge or face will always agree on the orientation of 
that edge or face. In addition to facilitating the treatment of tangential and normal 
traces, a unique orientation of edges and faces simplifies assembly of higher order 
Lagrange elements. A similar numbering scheme is proposed in the monograph [28j 
for tetrahedra in connection with //(curl) finite elements. Also, we note that the 
numbering scheme and the consistent facet orientation that follows render only one 
reference element necessary, in contrast to the approach of [2J. 
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Fig. 5.2. Two adjacent tetrahedra will always agree on the direction of a common edge tangent 
or face normal. The two tetrahedra in the figure share a common face defined by the global vertices 
V20, U50, and 1175. These three vertices have different local indices for K (2,3,4) and K' (1,2, 3), 
but the ordering convention, local numbering according to ascending global indices, ensures that 
both triangles agree on the direction of the common edges. In particular, the two tetrahedra agree 
on the directions of the first two edges of the common face and the direction of the face normal 
n oc ei X e2 = e\ X e' 2 . 

5.2. Mapping nodal basis functions. Next, we show how this numbering 
scheme and the FIAT choice of degrees of freedom give the necessary iJ(div) or 
_ff(curl) continuity. Assume that we have defined a set of nodal basis functions on 
K , that is, {$i}™ =1 such that 

£i(^j) = Sij, i,j = 1,2, ...,n, 

for a set of degrees of freedom These basis functions are mapped to two 

physical elements K and K' by an appropriate transformation T (contravariant or 
covariant Piola), giving a set of functions on K and K' , respectively. We demonstrate 
below that as a consequence of the above numbering scheme, these functions will 
indeed be the restrictions to K and K' of an appropriate global nodal basis. 

Consider H (curl) and a global degree of freedom t defined as the tangential com- 
ponent at a point a; on a global edge e with tangent t, weighted by the length of the 
edge e, 

£(v) — ||e|| v(x) ■ t = v{x) • e. 

Let ^ 7Curl be the covariant Piola mapping as before and let <p K and <p K be two basis 
functions on K and K' obtained as the mappings of two nodal basis functions, say, $ 
and on K , 

(t> K = J£ ,rl ($) and <j> K ' = ^ u , rl ($')- 
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Assume further that $ is the nodal basis function corresponding to evaluation of the 
tangential component at the point X € K along the edge E, and that $' is the nodal 
basis function corresponding to evaluation of the tangential component at the point 
X' e K along the edge E'. Then, if x = F K (X) = F K >{X'), the covariant Piola 
mapping ensures that 



Continuity for H(div) may be demonstrated similarly. 

In general, FFC allows elements for which the nodal basis on the reference element 
Kq is mapped exactly to the nodal basis for each element K under some mapping 
T, whether this be afHne change of coordinates or one of the Piola transformations. 
While this enables a considerable range of elements, as considered in this paper, it 
leaves out many other elements of interest. As an example, the Hermite triangle or 
tetrahedron [13] does not transform equivalently. The Hermite triangle has degrees 
of freedom which are point values at the vertices and the barycenter, and the partial 
derivatives at each vertex. Mapping the basis function associated with a vertex point 
value affinely yields the correct basis function for K, but not for the derivative basis 
functions. A simple calculation shows that a function with unit x-derivative and 
vanishing y-derivative at a point generally maps to a function for which this is not 
the case. In fact, the function value basis functions transform affinely, but the pairs 
of derivative basis functions at each vertex must be transformed together; that is, a 
linear combination of their image yields the correct basis functions. 

Examples of other elements requiring more general types of mappings include the 
scalar- valued Argyris and Morley elements as well as the Arnold- Winther symmetric 
elasticity element [7J and the Mardal-Tai- Winther element for Darcy-Stokes flow [37] . 
Recently, a special-purpose mapping for the Argyris element has been developed by 
Dominguez and Sayas 114 , and we are generalizing this work as an extension of the 
FIAT project as outlined below. 

If is the reference finite element basis and {cpf }, is the physical finite element 
basis, then equivalent elements satisfy <pf — Tk for each i. If the elements are 
not equivalent under Tk, then {</>f }i and {TK($i)}i form two different bases for 
the polynomial space. Consequently, there exists a matrix M K such that <pf = 
J2j Tk($j)- In the future, we hope to extend FIAT to construct this matrix M 
and FFC to make use of it in constructing variational forms, further extending the 
range of elements available to users. 

5.3. A note about directions. An alternative orientation of shared facets gives 
rise to a special case of such transformations. It is customary to direct edges in a 
fashion that gives a consistent orientation of the boundary of each triangle. However, 
this would mean that two adjacent triangles may disagree on the direction of their 
common edge. In this setting, normals would naturally be directed outward from 
each triangle, which again would imply that two adjacent triangles disagree on the 
direction of the normal on a common edge. It can be demonstrated that it is then 
more appropriate to define the contravariant Piola mapping in the following slightly 
modified form: 



<j) K {x) ■ e = $(A) • E = 1 and 
Thus, since e = e', it follows that 

t{<t> K ) = 



cj> K '(x')-e' = &(X')-E' = 1. 




T div {<S>) 



1 



Jk^oF- 1 ; 



|det J K \ 
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that is, the determinant of the Jacobian appears without a sign. 

To ensure global continuity, one would then need to introduce appropriate sign 
changes for the mapped basis functions. For two corresponding basis functions <\> K 
and <f) K as above, one would change the sign of <j) K or cf> K such that both basis 
functions correspond to the same global degree of freedom. Thus, one may consider 
obtaining the basis functions on the physical element by first mapping the nodal basis 
functions from the reference element and then correcting those basis functions with a 
change of sign: 

<f> K = ±t K . 

This would correspond to a diagonal M transformation where the entries are all ±1. 

Since a multilinear form is linear in each of its arguments, this approach cor- 
responds to first computing a tentative element tensor A K and then obtaining A K 
from A K by a series of rank one transforms. However, this procedure is unnecessary 
if the contravariant Piola mapping is defined as in (|2.9p and the numbering scheme 
described in section f5. II is employed. 

For nonsimplicial meshes, such as meshes consisting of quadrilaterals or hexa- 
hedra, the situation is somewhat more complicated. It is not clear how to ensure 
a consistent, common local and global direction for the edges. Therefore, the UFC 
specification instead requires a consistent orientation of the boundary of each cell. In 
this situation, the alternative approach, relying on the introduction of sign changes, 
is more appropriate. 

6. Examples. In order to demonstrate the veracity of the implementation and 
the ease with which the -ff (div) and -ff (curl) conforming elements can be employed, 
we now present a set of numerical examples and include the FFC code used to de- 
fine the variational forms. In particular, we return to the examples introduced in 
section [21 which include formulations of the Hodge Laplace equations, the cavity res- 
onator eigenvalue problem, and the weak symmetry formulation for linear elasticity. 

6.1. The Hodge Laplacian. Consider the weak formulations of the Hodge 
Laplace equation introduced in Examples 12.11 and 12.21 For Q, C R 2 and differential 1- 
and 2-forms, we have the mixed Poisson equation (I2.ip . Stable choices of conforming 
finite element spaces S/, X Vf, C iJ(div) x L 2 include Vh = DG r _i in combination 
with Y> h e {RT r _i, BDFM r , BDM r } for r = 1, 2, . . . . The FFC code corresponding 
to the latter choice of elements is given in Tabic [B~T1 Further, for ft C R 3 , we give the 
FFC code for the formulation of (|2 . 3[) with the element spaces NED*_i x RT r _i C 
ff(curl) x ff(div) in Table O 

For testing purposes, we consider a regular tessellation of the unit square/cube, 
fi = [0, l] n , n = 2,3, and a given smooth source for the two formulations. In particu- 
lar, for (|2.1[) . we solve for 

(6.1) u(x\,X2) = Csin(7rcci) sin(7TX2), 

with C a suitable scaling factor, and for (|2.3p . we let 

(x1(xi — l) 2 sin(7TX2) sin(7ra;3)\ 
cc2(a;2 - l) 2 sin(7rxi) sin(7ra;3) . 
xf(xs — I) 2 sin(7rxi) sin(7riE2) / 
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Table 6.1 

FFC code for the mixed Poisson equation. 



r = 3 

S = FiniteElementO'BDM" , "triangle", r) 
V = FiniteElementO'DG" , "triangle", r - 1) 
element = S + V 

(tau, v) = TestFunctions(element) 
(sigma, u) = TrialFunctions (element) 

a = (dot (tau, sigma) - dot(div(tau) , u) + dot(v, div (sigma) ) *dx 
L = dot(v, f)*dx 



Table 6.2 

FFC code for the curl-div formulation of the Hodge Laplace equation. 



r = 2 

CURL = FiniteElementC'Nedelec" , "tetrahedron", r - 1) 
DIV = FiniteElementC'RT" , "tetrahedron", r - 1) 
element = CURL + DIV 

(tau, v) = TestFunctions(element) 
(sigma, u) = TrialFunctions (element) 

a = (dot (tau, sigma) - dot (curl (tau) , u) + dot(v, curl (sigma)) \ 

+ dot(div(v), div(u)))*dx 
L = dot(v, f)*dx 



Note that u given by (|6.2j) is divergence-free and such that u x n = on the exterior 
boundary, and thus satisfies the implicit natural boundary conditions of (|2.3j) . 

A comparison of the exact and the approximate solutions for a set of uniformly 
refined meshes gives convergence rates in perfect agreement with the theoretical values 
indicated by Table I3.2[ up to a precision limit. Logarithmic plots of the L 2 error of 
the flux using E/j € {RTV_i, BDM r } versus the mesh size for r = 1, 2, . . . , 7 can be 
inspected in Figure HO for the mixed Poisson problem (with C = 100). 

For the curl-div formulation of the Hodge Laplace equation (|2.3[) . we have included 
convergence rates for u and a in Table 16.31 Note that the convergence rates for the 
combinations NED,, x RT r and NED r x BDM r , r = 1, 2, are of the same order, except 
for the || • ||div error of it, though the former combination is computationally more 
expensive. 

6.2. The cavity resonator. The analytical nonzero eigenvalues of the Maxwell 
eigenvalue problem (|2.5j) with O = [0, ir] n , n ~ 2, 3, are given by 

(6.3) uj 2 = m \+ml-\ Vm 2 nl vrn G {0} U N, 

where at least n — 1 of the terms mi must be nonzero. It is well known [28^ that 
discretizations of this eigenvalue problem using H 1 conforming finite elements produce 
spurious and highly mesh-dependent eigenvalues lu 2 . The edge elements of the Nedelec 
type, however, give convergent approximations of the eigenvalues. This phenomenon 
is illustrated in Figure 16.21 There, the first 20 nonzero eigenvalues (J 2 . N produced 
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Fig. 6.1. Convergence rates for the discretized mixed Poisson equation i2.lt using RT r _i X 
DGr-i (left) and BDM r X DG r —i (right), r = 1,2, ... ,7. Logarithmic plots of the L 2 error of the 
flux approximation: \\a — (T^||o versus mesh size. The convergence rates in the left plot are 0(h r ) 
and the convergence rates in the right plot are 0(h r+1 ); cf. Taftle EOl The error does not converge 
below ~ 10~ 10 in our experiments as a result of limited precision in the evaluation of integrals 
and/or linear solvers. The exact source of the limited precision has not been investigated in detail. 



Table 6.3 

Averaged convergence rates for the discretized curl-div formulation of the Hodge Laplace equa- 
tion IpO) using NED r _i X RT r _i, r = 1, 2, 3, and NED r X BDM r , r = 1, 2. Number of degrees of 
freedom in the range 80,000-300,000. 



Element 


k - u h Ho 


k ~ a h curl 


II" - "hllo 


Ik - "felldiv 


NED x RT 


0.99 


0.98 


0.99 


0.98 


NEDi x BDMi 


1.96 


2.00 


1.95 


0.96 


NEDi x RTi 


1.97 


1.97 


1.98 


1.98 


NED 2 x BDM 2 


3.00 


2.99 


2.97 


1.97 


NED 2 x RT 2 


2.98 


2.96 


2.97 


2.97 



18- 

17 

16 



• • Nedelec 
a a Lagrange 



10 - 
9 

8 - 



Fig. 6.2. The first 20 eigenvalues of the cavity resonator problem computed using first order 
Nedelec elements (NET)q) and Lagrange elements (Pi) on a coarse (16 X 16) criss-cross mesh. The 
exact analytical values are indicated by the horizontal grid lines. 
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Table 6.4 

FFC code for linear elasticity with weak symmetry. 



def A(sigma, tau, nu, zeta) : 

return (nu*dot ( sigma, tau) - zeta*trace (sigma) *trace (tau) ) *dx 

def b(tau, w, eta): 

return (div(tau[0] )*w[0] + div(tau[l] )*w[l] + skew (tau) *eta) *dx 

nu = 0.5 
zeta = 0.2475 
r = 2 

S = FiniteElementO'BDM" , "triangle", r) 

V = VectorElementC'Discontinuous Lagrange", "triangle", r-1) 
Q = FiniteElementC'Discontinuous Lagrange", "triangle", r-1) 
MX = MixedElement ( [S , S, V, Q] ) 

(tauO, taul, v, eta) = TestFunctions (MX) 
(sigmaO, sigmal, u, gamma) = TrialFunctions(MX) 
sigma = [sigmaO, sigmal] 
tau = [tauO, taul] 

a = A(sigma, tau, nu, zeta) + b(tau, u, gamma) + b(sigma, v, eta) 
L = dot(v, f)*dx 



by the Nedelec edge elements on a regular criss-cross triangulation are given in com- 
parison with the corresponding Lagrange eigenvalue approximations u\ L . Note the 
treacherous spurious Lagrange approximations such as ui\ L « 6, 15. 

6.3. Elasticity with weakly imposed symmetry. As a final example, we 
consider a mixed finite element formulation of the equations of linear elasticity with 
the symmetry of the stress tensor imposed weakly as given in Example 12.41 In the 
homogeneous, isotropic case, the inner product induced by the compliance tensor A 
reduces to 



for i/, £ material parameters. A stable family of finite element spaces for the discretiza- 
tion of ([H]) is given by [6 : BDM 2 x BG 2 r _ 1 x DG r _i C iJ(div; fi,M) x L 2 (fi,R") x 
L 2 (f2), r = 1, 2, . . . . The 17 lines of FFC code sufficient to define this discretization 
are included in Table 16.41 

Again to demonstrate convergence, we consider a regular triangulation of the unit 
square and solve for the smooth solution 



The theoretically predicted convergence rate of the discretization introduced above is 
of the order 0(h r ) for all computed quantities. The numerical experiments corrobo- 
rate this prediction. In particular, the convergence of the stress approximation in the 
_ff(div) norm can be examined in Figure [6~31 



(r,Aa) 



v{t, a) — C(tr t, tr a) 



(6.4) 
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Fig. 6.3. Left: Convergence rates for elastic stress approximations of l|2.6|l . Logarithmic plot 
of H (div) error of the approximated stress a versus mesh size. The convergence rates are 0(h r ), 
r = 1,2,3,4. Right: Elastic dolphin hanging by the tail under a gravitational force. 

7. Conclusions. The relative scarcity of if (div) and -ff (curl) mixed finite ele- 
ment formulations in practical use may be attributed to their higher theoretical and 
implcmentational threshold. Indeed, more care is required to implement their finite 
element basis functions than the standard Lagrange bases, and assembly poses ad- 
ditional difficulties. However, as demonstrated in this work, the implementation of 
mixed finite element formulations over iJ(div) and iJ(curl) may be automated and 
thus be used with the same ease as standard formulations over H , In particular, the 
additional challenges in the assembly can be viewed as not essentially different from 
those encountered when assembling higher order Lagrange elements. 

The efficiency of the approach has been further investigated by 01gaard and 
Wells [32], with particular emphasis on the performance when applied to more com- 
plicated PDEs. They conclude that the tensor representation significantly improves 
performance for forms below a certain complexity level, corroborating the previous re- 
sults of [21]. However, an automated, optimized quadrature approach, also supported 
by FFC, may prove significantly better for more complex forms. These findings indi- 
cate that a system for automatically detecting the better approach may be valuable. 

The tools (FFC, FIAT, DOLFIN) used to compute the results presented here are 
freely available as part of the FEniCS project [15] and it is our hope that this may 
contribute to further the use of mixed formulations in applications. 
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